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Abstract. We are concerned with the numerical simulation of interaction of 
incompressible flow with elastic bodies. Our goal is to simulate airflow in human 
vocal folds and their flow-induced vibrations. We consider two-dimensional viscous 
incompressible flow in a time-dependent domain. The fluid flow is described by the 
Navier-Stokes equations in the arbitrary Lagrangian-Eulerian formulation. The flow 
problem is coupled with the elastic behaviour of the solid bodies. The dynamical 
problem of deformation of the elastic body is formulated. The developed solution 
of the coupled problem based on the finite element method is demonstrated by 
numerical experiments. 


Introduction 


In many cases we solve either problems of the fluid mechanics or problems of the solid 
mechanics. But in many applications we are interested in the behaviour of both the fluid flow 
and the solid and moreover their interaction. A lot of examples we find in biomechanics. Here 
we are concerned with the phonation onset in biomechanics of voice, which is an important 
characteristic of human voice production. The phonation onset can be characterized as a state 
of the system when it is loosing the aeroelastic stability, i.e. when the airflow parameters 
like subglottal pressure or airflow rate cross a limit value and the system becomes unstable 
by flutter. Frequency-modal analysis of a simplified three-mass model of the vocal fold in 
interaction with a potential flow separated at the superior edge of the vocal fold showed that 
the two eigenfrequencies are coupled when the instability occurs [Horáček et al., 2002]. A similar 
linear stability analysis was performed on a 2-D continuum model of vocal folds in potential 
flow by [Zhang, 2009]. Here the instability threshold is studied in the time domain using 2-D 
FE model of the vocal folds coupled with 2-D FE model of viscous incompressible flow. 


Fluid flow in a moving domain 


We consider incompressible viscous flow in a bounded domain of C R? depending on time 
t € [0,7]. By v = v(a,t) we denote the velocity and by p = p(a,t) the kinematic pressure (i.e., 
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Figure 1. The model of vocal fold. 
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pressure divided by the density o” of the fluid), x € al , t € (0,7) and v denotes the kinematic 
viscosity. Incompressible viscous flow is described by the system of the Navier-Stokes equations 
equipped with initial and boundary conditions [Feistauer, 1993]. 


—+(v-V)vu+Vp—vAv=0 inff, (1) 
V-v=0 ino, (2) 
In order to simulate flow in a moving domain, we employ the Arbitrary Lagrangian-Eulerian 


(ALE) method. This method is based on a special mapping of the reference configuration of 


onto the deformed, actual configuration al . We reformulate the Navier-Stokes equations in the 
ALE form [Nomura, 1992]: 


—v4+((v—w):-V)v+Vp—vAv=0 in Qf, (3) 

V-v=0 inl, (4) 

Here pe is so-called ALE derivative and w denotes the domain velocity. System (3) - (4) is 
equipped with the initial condition 

v(z,0) =v, £E of, (5) 


and boundary conditions. We assume that anf = A U r U Tw,, where Ei ie and Ty, are 
mutually disjoint. On ph representing the inlet and impermeable fixed walls, we prescribe the 
Dirichlet boundary condition, on the impermeable moving walls [w+ we assume that the fluid 
velocity v is equal to the domain velocity of the elastic body and on the outlet ine we prescribe 
the so-called ”do-nothing” boundary condition: 


Ups = vp oI, (6) 
rw = wiry, olw,, (7) 
Ov 
=(P- Pref) + Uz = 0 on rÊ. (8) 


Here n is the unit outer normal to onf and pref is a prescribed reference outlet kinematic 
pressure. 

There are several possibilities how to carry out the space-time discretization. In order to 
develop a stable, accurate scheme, which can easily treat complicated boundaries, we apply the 
finite element method (FEM). Here, the Taylor-Hood P2/P, finite element pair satisfying the 
Babuska-Brezzi condition is used. In order to avoid spurious oscillations in approximate solu- 
tions in the case of high Reynolds numbers we apply the streamline diffusion method together 
with div-div stabilization of the pressure. For the time discretization we use a second-order 
two-step backward difference formula using the computed approximate solution v”~! in al 


n—1 
and v” in of for the calculation of v”*t in the domain of ,,- Here the solutions from pre- 
vious time steps are transformed via ALE mapping on the actual computational domain. For 
the numerical solution of incompressible flow we use the program FEMFLUID developed by P. 
Sváček. The space-time discretization is specified in [Sudéek, 2004]. 


The linear problem of elasticity 


In what follows, Q? c R? will be a bounded domain representing the elastic body of the 
form of vocal folds. We denote by u(a,t), x € 9°, t € (0,7), the displacement of the body and 
the strain tensor as 


a 1 Ou; Ou; ay Neo 
est) = 5 (grt ge), 65-12, (9 


62 


KOSIK: NUMERICAL SIMULATION OF INTERACTION OF FLUID FLOW AND ELASTIC BODY 


The deformation of the vocal folds is modelled by the generalized Hooke’s law for isotropic 
bodies [Nečas et al., 1981] 


T = Adiv Udi; = 2Ueij, i,j = 1,2, (10) 


where GAH j=1 denotes the stress tensor and A and p are the Lamé coefficients related to the 


Young modulus E and to the Poisson ratio o as 


p(3A + 2u) À 
E = => — 11 
Naa ETT an 


The dynamic equations of an elastic body have the form 


3u; a 
Fs pOUi b 
op +Co À — 0, on (0,T)x Q’, i=1,2. (12) 


b oui 


The expression Co’, where C > 0, is a dissipative damping of the system and o? denotes 
the density of the solid nigeria: We complete the elasticity problem with initial and boundary 
conditions. The initial conditions read 


u(-,0)=0, =—(,0)=0, mm. (13) 


Further, let Q? = Tw U gee where I'w and ry are two disjoints parts of 00°. Let surface force 
T” be prescribed on the boundary Tw and let the part of the boundary T$ be fixed: 


Sty SP only xk OP), 2=12, (14) 


u = 0 onI, x (0,7). (15) 


We are looking for the displacement u satisfying equation (12) and the initial and boundary 
conditions (13) - (15). 

Further, we carry out the space-time discretization. We reformulate the problem in a weak 
sense. For given u? € H!(°), z° € L?(0°) we seek for the solution u € L? (0, T; V), if it holds 
u' € L? (0,T; £7(0°)), u” € L? (0,T; V*), 


q? d 5 
qe Y u(t), Yoo» += J (Co’ult), y)o,o +alu, y;t) = (T” (t) Worw, Yy €V, te [0,7], 


V=V?, Cae) | lr, = 0}, 


OYi 
alu, y;t) =f, Sí (Adivu(t)d;; + 2ue;; (u )) 55, dæ. (16) 


i, j=1 
We apply the finite element method using continuous piecewise linear elements. Thus we seek 
the approximate solution up on the space L? (0, T; V;,), where V, = VŽ, 


V, = fon c C(%) | ọ is linear on each triangle of the triangulation, yrs, = o} CV. (17) 


The semi-discretized problem can be written as a second order system of ordinary differential 
equations. For the time discretization we apply the Newmark scheme. In each time-step we get 
a linear algebraic system with symmetric positive definite matrix. The solution of this system 
was realized by the solver, which is based on the conjugate gradient method. 
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The coupled problem 


Up to now we assumed that the fluid flow and the deformation of the elastic body were 
separated processes. On the common boundary of both domains we need to take into account 
mutual action of the fluid and the body. We denote as the common boundary set [w,, defined 
as 

Tw, = {x € R°; x = X +u(X,t), X eTy}. (18) 


The domain Qf is defined by the displacement u of the part Tw at time t. The ALE mapping 
A; is defined with the aid of a special stationary linear elasticity problem. If we know the 
computational domain al obtained by the fluid at time t, we can solve the problem describing 
the flow and we can assign the surface force T” acting on the body 2° on the part of the 
boundary Tw. This surface force is opposite to the force Tf acting by the airflow on the 
reference domain Of, 

T” = -T'. (19) 


We obtain the force T/ by the transformation of the force Ti defined on D by the equation 
2 
TE eye, (20) 
j=1 


where 7 is the unit outer normal to the boundary ['w,, and te are components of the stress 


tensor Tf. We can compute the stress in the actual configuration Oh by the equation: 


ðv; Ov; me 
i =e Gaz ast (= F 2) ERTE ie 
J a 


The ALE mapping describing the movement of the computational domain matches to the def- 
inition of deformation. Thus we use Piola’s transformation for the acting force TÍ onto the 
boundary of the reference domain of : 


TS = |cof (VA, )n| TS (22) 


where cof(V.A;,,) denotes the cofactor matrix of VA, . Since n = St we get 


TS = rf cof (VAs, )n. (23) 


The procedure of the implementation is following. First we approximate the computational 
domain of by the domain of i- The computation of the flow problem is carried out in the 


approximate domain oe We get the approximate velocity vp and pressure pp. We use the 


obtained result to compute an approximation of the components of stress tensor re on the 


I'y,,,- The stress tensor gives us the surface force Ti on the boundary Ty,,. We transform 
the surface force TÍ to obtain the surface force Tf acting on the reference domain Q? on the 
part of boundary I'w. Further we solve the problem of elasticity and get the approximate 
displacement up at time tn. Hence we are able to solve the flow problem in the updated domain 
‘on and repeat the computation of the elasticity problem with the updated approximation of the 
stress tensor Tey. If the approximation of the displacement up, changes more than a prescribed 
tolerance, we repeat the whole process to obtain a better approximation of the flow problem 
solution and consequently a better approximation of the body displacement. Then we come to 
the next time step tn+1- 
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Figure 2. Material characteristics. 


Numerical experiments 


We consider model of human vocal folds and vocal tract [Feistawer et al., 2010]. The 
vocal folds represented by the domain Q? have different material characteristics in different 
subdomains. We use the same time step T = 5-107°s for the solution of the coupled flow and 
elasticity problem and the input data v = 1.5 - 1075 m?.s~!, of = 1.17 kg.m™’, C = 0.1 s7}, 
the initial velocity v? = 0 and the inlet and outlet pressure pin = 600 Pa, Pout = 0 Pa. 

As for the domain Q? we distinguish subdomains with different material characteristics. 
We assume that the material density o? = 1040 kg.m~° is the same for all subdomains, but the 
values of the Young modulus E and the Poisson ratio o are different, see Table 2 and Figure 1. 
On the part of boundary I'w we prescribe the conditions for the moving boundary. We acquire 
the surface forces T” by solving the problem of fluid flow. The computational process starts by 
the solution of the flow problem in the domain Qf at the initial time tg = —2.1073 s. At time 
t = 0 the structure was released and the solution of the complete fluid-structure interaction 
started. 

Figure 3 shows velocity streamlines and displacement of the computational domain at 
several time instants. 


Conclusion 


The numerical method for solving the flow-induced vibrations of an elastic body in incom- 
pressible viscous flow has been developed and applied to the simulation of airflow in interaction 
with human vocal folds. The results are in good agreement with previous simulations using 
simplified models and with known physiological data (see, e.g. |Hordéek, 2002]). 

In future we will focus on the verification of our computations. We are going to perform 
more numerical experiments with various time steps and various number of finite elements. 
Another goal of our future work will be the simulation of the complete closure of the channel 
and the use of a non-linear elasticity model in the coupled problem. 


Acknowledgments. This research was supported by the grant SVV-2011-263316 of the Charles 
University in Prague. 


References 


Feistauer, M.: Mathematical Methods in Fluid Dynamics, Longman Scientific & Technical, Harlow, 1993. 


65 


KOSIK: NUMERICAL SIMULATION OF INTERACTION OF FLUID FLOW AND ELASTIC BODY 


Feistauer, M., Cesenek, J., Horáček, J., Kučera, V. & Prokopová, J., DGFEM for the numerical solution 
of compressible flow in time dependent domains and applications to fluid-structure interaction. In: 
ECCOMAS CFD 2010, J.C.F. Pereira and A.Sequeira (Eds), Lisbon, Portugal, CDROM, ISBN 978- 
989-96778-1-4, 14-17 June 2010. 

Horáček, J., Švec, J. G.: Aeroelastic model of vocal-fold-shaped vibrating element for studying the 
phonation, Journal of Fluids and Structures 16 (7), 931-955. doi:10.1006/jfls.454, 2002 

Nečas, J. & Hlaváček, I., Mathematical Theory of Elastic and Elasto-Plastic Bodies, An Introduction. 
Elsevier, Amsterdam, 1981. 

Nomura, T. & Hughes, T. J. R.: An arbitrary Lagrangian-Eulerian finite element method for interaction 
of fluid and a rigid body, Comp. Methods Appl. Mech. Engrg., no. 95., 115-138, 1992. 

Sváček P. & Feistauer M.: Application of a stabilized fem to problems of aeroelasticity. In M. Feis- 
tauer, V. Dolejsi, and Najzar K., editors, Numerical Mathematics and Advanced Applications, ENU- 
MATH2003, pages 796-805, Heidelberg, Springer, 2004. 

Zhang, Z.: Characteristics of phonation onset in a two-layer vocal fold model, J. Acoust. Soc. Am. 
125(2), 1091-1102, 2009. 


Figure 3. Streamlines at time instants t = 0.033, 0.034, 0.036 and 0.03725 s 
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